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ABSTRACT 

An algebraic adaptive grid system based on equidistribution law and utilized by the Non-Uni- 
form Rational B-Spline (NURBS) surface for re-distribution is presented. A weight function, utiliz- 
ing a properly weighted boolean sum of various flow field characteristics is developed. Computa- 
tional examples are presented to demonstrate the success of this technique. 

INTRODUCTION 


Grid adaptive schemes are divided into two basic categories: differential and algebraic. The dif- 
ferential method is based on a variational approach where a function which contains a measure of 
grid smoothness, orthogonality, and volume variation is minimized by using a variational principle. 
This approach provides a solid mathematical basis for the adaptive method, but the Euler-La- 
grange equations must be solved in addition to the original governing equations. On the other hand, 
the algebraic method requires much less computational effort, but the grid may not be smooth. The 
algebraic techniques are based on devising an algorithm where the grid movement is governing by 
estimates of the local error in the numerical solution. This is achieved by requiring the points in the 
large error regions to attract other points and points in the low error region to repel other points. One 
of the disadvantages of algebraic method is that it generates the unsmooth grid. The utilization of 
NURBS will generate the well-distributed smooth grid (Fig. 1). 

An adaptive grid is a grid that controls the placement of grid points automatically based on the 
solution of the physical problem under consideration and allows optimal grid redistribution as the 
solution progresses. The accuracy of the numerical algorithm depends not only on the formal order 
of approximation but also on the distribution of grid points in the computational domain. The adap- 
tive grid is one of the methods that can make the numerical algorithm more accurate. The grid adap- 
tive scheme, presented in this paper, is based on the equidistribution principle. 
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The development of the adaptive algorithm for the structured flow simulation is accomplished 
as a two step process. The first step is to define an adaptive weighting mesh (distribution mesh) (Ref. 
1 -2) on the basis of the equidistribution law applied to the flow field solution. The second, and prob- 
ably the most crucial step, is to redistribute grid points in the computational domain according to 
the aforementioned weighting mesh. Adaptive weighting (distribution mesh) provides the informa- 
tion on the desired concentration of points to the grid redistribution scheme. The evaluation of 
weighting mesh is accomplished by utilizing the weight function representing the solution vanation 
and the equidistribution law (Ref. 3). The selection of the weight function plays a key role in gnd 
adaptation (Ref. 4-5). A new weight function utilizing a properly weighted boolean sum of various 
flowfield characteristics is defined. The redistribution scheme is developed by utilizing Non-Uni- 
form Rational B-Spline Surface (NURBS) representation (Ref 6). The application of NURBS rep- 
resentation results in a well-distributed smooth grid by maintamg the fidelity of the geometry 
associated with boundary curves. Various computational examples are presented to demonstrate the 
success of these methods. 


GRID ADAPTATION METHOD 


The basis of the grid adaptation is the equidistribution law in one dimension (Ref. 8) applied to 
the flow field. The equidistribution law is represented as 

x t w = cons Unit ( 1 ) 

With the weight function w taken as a function of £, this is just the Euler equation for the minimi- 
zation of the integral 

/ = ( w ( n x* <% (2) 

Applying the spring analogy form to the above integral, the point location x, can be defined by 
(Ref. 8) 


x i 



( i— 2, 3 , N— 1 ) 


(3) 


The integral in the denominator depends on the point distribution, amounting to a sum of 1/w 
over the points. However = 1 by construction regardless of the distribution. Therefore, the inte- 
gral in (3) is evaluated iteratively. The integral in the numerator is then also re-evaluated for each 
point, thus changing the point distribution again. This process must be continued until convergence 
before the final new distribution is obtained. 
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Figure 1. 

Use NURBS surface to generate smooth grid. 




ALGORITHM 


A self-explanatory, pictorial view of the algorithm of algebraic adaptive gnd generation is 
shown in Fig. 1 . 

The grid and solution will be transferred from the flow solver if the solution is good enough to 
apply grid adaptation. The criterion of good solution can be the pressure L 2 residue or another prop- 
erty of solution that can be used as the convergent indication. 
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Figure 2. 

Algorithm for 2-D Algebraic Adaptive Grid Generation System 








WEIGHT FUNCTION 


The weight function is a very important part in the adaptive grid system. Here we introduce a 
new weight function applicable to various flow field characteristics. The weight functions are com- 
puted in both the ^ and r\ computational directions, and then coupled adaptation is applied. A linear 
combination 


1 + Xj Wj , with y* Aj = 1 (4) 

where N - number of flow parameter ( e.g. pressure, temperature, density, etc.) 

Xj - weighting factor associated with flow parameter Xj 2 0 

wj - ctj qj © Pj kj = ctj qj + pj kj - («j + pj -1) qj kj 

qj - scaled gradient of the flow variable j such that 0 s qj s 1 

kj - scaled curvature values of the flow variable j such that 0 s kj ^ 1 

0 s ctj s 1, 0 ^ pj s 1,- weight factors assigned to gradients and curvatures is developed as 
the weight function utilizing the Boolean sum of contributions from scaled gradients and curvatures. 
The Boolean sum allows an appropriate weight to the gradients as well as the curvatures. For exam- 
ple, in the case of the 7 th flow parameter, the influence of weights can be presented as follows: 
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Note that the value of the weight contribution is controlled by the weight factors and is at a maxi- 
mum when gradients and/or curvature values are at a maximum. The gradients and curvatures 
associated with the flow characteristics under consideration are smoothed using a four point diffu- 
sion process. 
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SCALING SCHEME FOR a & p 


Dwyer has developed a method to determine the percentage change in a dependent variable as 
a priori. ( Ref. 9 ). This method is enhanced here to determine the suitable scaling a and P at each 
grid point. (Ref. 3) The adaptation technique can be described in the generalized coordinates as 


£ = 


ai (]JZ+ fil fair - (a + /* - 1 ) f tjiiff 

Jp JO Jo Jo 

a f qdx + ft f kdx — (c* + ft — 1) f qdx f 

Jo Jo Jo Jo 


(5) 


kdx 


where x - normalized arc length, q - flow parameter gradient, k - flow parameter curvature, a & 
p - weight factor assigned to gradients and curvatures, so that 


M = 


[ aq + ftk — (a + ft - 1 )qk ] Ax 


a I qdx + ft I 

Jo Jo 

Then with Rj defined as 


kdx - (a +0 - 1) 
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we have the maximum percentage change in the solution over a grid interval. 
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since AS; = 1/N and N+l is the number of points on the coordinate line. If we take an equality in 
Eq. (8) with Ri from Eq. (7) , we have 


NR, 
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Ignoring the effect of the p term: 
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Similarly with R 2 defined as 


*2 = 
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we can get 
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(13) 


With Boolean sum operation, we found that the choice of Rj = 1 and R 2 = 1 still results in a very 
good adaptive grid. 


GENERATE REDISTRIBUTED ALGEBRAIC GRID 


All of the adaptive grids displayed in this paper are generated by the Non-Uniform Rational B- 
spline in surface (NURBS) or reparametrized uniform B-spline surface. The reparametrized uni- 
form B-spline curve is applied to move the boundary points in order to keep the sharp corners 
(Fig.3). In general, the degree m NURBS surface (3-D) is defined as the projection of a tensor prod- 
uct B-spline in 4-D (Ref. 5-6): 


r <$,*) 




S S <M) h , 

»»0 ; =0 


(14) 


where are the mth order bivariate B-spline basis functions. H,j and Q,j are the associated latice 
of weights and control points. In the grid redistribution process, a matrix form of bi-cubic NURBS 
was applied which can be formulated as (Ref. 7): 
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where control points and weights [ H Q ] are as follow: 
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The matrix M h is defined as 
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The redistributed algebraic grid is generated by utilizing distribution mesh as parameter space. 
The convex hull, local support, and variation diminishing properties (Ref. 7) of B -spline functions 
contribute to the generation of the well-distributed smooth grid. The application of inverse NURBS 
formulation (Ref. 6) allows reevaluation of control points which influences the fidelity of surface 
geometry during the redistribution process. The redistributed grid is evaluated for positive areas and 
smoothness. If there is any negative area, the further more grid adaptation will not be executed. 
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APPLICATIONS 


The initial gnds were generated by (GENIE-2D-3D) (Ref. 1). The basic flow simulanon codes 
is the computer code PARC2D (Ref. 10) which are algorithmically based on the Beam-Warming 
implicit finite difference scheme using approximate factorization. The routines for gnd adaptation 
are incorporated in the PARC2D system with proper data flow. 

( 1) Wedge Problem with a 5 Degree Incline Angle at the Leading Edge. Using the theory devel- 
oped for supersonic flow, it is possible to compute the exact behavior of the shock waves generated 
by single and double wedge in supersonic flow. In the wedge problem, the Mach number is 1 .9, and 
the Reynolds number is 1.0 x 10 4 * * * * * 10 with a grid size at 80 x 60. There is one shock and one reflected 
shock in the control volume. Without grid adaptation, the shock and the reflected shock are as thick 
as two small tapes. Even with a significant number of iterations (e.g. 10000 iterations)(Fig. 4). of the 
flow solver, the solution shows very little improvement. After grid adaptation, the two shocks appear 
thinner. Extremely thin shocks can be achieved by increasing grid adaptation numbers (Fig. 5). The 
grid adaptation was performed at every 300 iterations. The redistributed gnds and their influence 
on the solution are presented in Figure 5. It can be seen that at every adaptation stage the pressure 
contours have indicated a crisper shock. A closer look at the changing gnd behavior near the shock 
region is presented in Figure 6-7. In spite of extremely skewed grid lines, the code was able to 
execute and provide a better solution. 

(2) Double Wedge Case with Grid Size 103 x 27 and 2.0 Mach Number. As expected, the shocks 
are thinner with grid adaptations, and the outlet near the mixture area of the second shock and the 
third shock is more clear (Fig. 8). 

(3) Cylinder Case with Grid Size 80 x 100, Mach Number 3.0, and Zero Angle of Attack. The 
bow shock in front of the cylinder is captured better in the case of the adapted grid as shown in figure 
9-10. Although the bow shock is a curve, the adaptive grid also bends along the shock wave. 

It is generally well known that either rapidly varying the step size or the skew grid angle can 
lead to poor accuracy when using finite difference approximations. Dwyer has shown by numerical 
experiments that variable grids coupled with adaptive techniques are accurate. (Ref. 9) In this paper 
we also observed in several cases that high skew angle grids coupled with adaptive techniques can 
excellently resolve the shock (Fig. 5-10). We also acknowledge that if the orientation of the gnd 
line is aligned with the shock wave, the shock wave is captured exactly (Fig. 11). Thus, not only the 
cluster of the grid points near the shock is necessary but also the alignment of the grid points with 
the shock wave is important. We found that during the adaptation, the alignment of the gnd points 
with shock wave will be accomplished by forming the skewed grid lines, according to the flow solu- 
tions. It is hoped that mathematical analysis could be developed to prove this observation, but the 
present investigators have not proceeded in this pursuit. We hope it will be accomplished in the near 
future. 

( 4) Airfoil NACA 0012 with Grid Size is 180 x 160, Mach Number 0.799, and the Angle of Attack 

1.25. The upper surface shock is thinner with the grid adaptation. The lower surface has no weak 

shock if no grid adaptation is applied, but with grid adaptation, the weak shock was captured (Fig. 

1 2) . Comparing the C p vs x/c charts, it is clear tharthe adaptive grid caught both shocks much better 

than the initial grid did (Fig. 13), but the adaptive grid loses the grid concentration at the leading 

edge and the tail. This is because the control field is so large that the effect of concentrating grid lines 

near the high gradient and/or the high curvature is spread to the whole field. The strategies that can 
improve this kind of condition are modifying the weight function by controlling the distribution near 
the body surface or subdividing the whole field into several blocks with some smaller blocks near 
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the surfaces of the airfoil and applying the local grid adaptation. The second strategy needs the devel- 
opment of a multi-block adaptive grid system. Both strategies are currently under investigation. 


CONCLUSION 


With a new weight function formula, a new method that applies the NURBS surface in the alge- 
braic adaptive grid generation system is introduced. The effect of this new adaptive grid system is 
pretty good in capturing the shock wave for some smaller control field. The grid adaptation strategies 
that are applied to catch the high gradient and/or high curvature area are used not only to cluster grid 
points near those areas but also to move the grid points to align along the high solution derivative 
areas. For complicated control fields, the multi-blocks adaptive grid system may be required. The 
weight function that is suitable for laminar or turbulance flow simulation which need very small cells 
near body surface, also needs to be investigated. The extension of this 2-D general adaptive grid 
system to a three dimensional problem and its adaptation to a multiblock environment is currently 
underway. 
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Figure 4. 

Single wedge initial grids & solutions after 300 & 10000 iterations. 
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Figure 5. 

Single wedge grid adaptation process. 
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Figure 6. 

Single wedge, closer pressure contour look near shock region. 





flrWtfd 



Figure 7. 

Single wedge, closer Mach number contour look near shock region. 







Figure 8. 

Double wedge case non-adaptive/adaptive simulation. 
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Figure 9. 

Flow around cylinder initial/adaptive cases. 
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Figure 11. 

Grid lines align with the shock wave capture shock exactly. 
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Figure 12. 

Airfoil NACA0012 case non-adaptive/adaptive grid/solution. 



Cp vs X/C on airfoil 0012 for Initial A Adaptive Solutions 

lfaeh *= 0.799, ife>ha b f.ZS, Hr * 9000, grid *i»« b 190x100 



Airfoil NACA0012 case Cp vs x/c plot. 






